#ifndef __BWLIB_NUMBER_THEORY_GCD_H__
#define __BWLIB_NUMBER_THEORY_GCD_H__
#include "../config.h"
#include <utility>
#ifdef __BWLIB_WITH_GMPXX
#include <gmpxx.h>
#endif

namespace bwlib
{
	/* @brief: 计算 gcd(a, b)
	   @complexity: O(log n)
	   @param: a, b
	   @return: gcd(a, b) */
	template<typename Type>
	inline Type gcd(const Type& a, const Type& b)
	{
		Type r, ta = a, tb = b;
		while(b > 0)
		{
			r = ta % tb;
			ta = __BWLIB_MOVE(tb);
			tb = __BWLIB_MOVE(r);
		}
		return __BWLIB_MOVE(ta);
	}

	/* @brief: 计算 lcm(a, b)
	   @complexity: O(log n)
	   @param: a, b
	   @return: lcm(a, b) 
	   @warning: integer overflow */
	template<typename Type>
	inline Type lcm(const Type& a, const Type& b)
	{
		return a / gcd(a, b) * b;
	}

#ifdef __BWLIB_WITH_GMPXX
	template<>
	inline mpz_class gcd(const mpz_class& a, const mpz_class& b)
	{
		mpz_class g;
		mpz_gcd(g.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
		return g;
	}

	template<>
	inline mpz_class lcm(const mpz_class& a, const mpz_class& b)
	{
		mpz_class g;
		mpz_lcm(g.get_mpz_t(), a.get_mpz_t(), b.get_mpz_t());
		return g;
	}
#endif

#ifdef __BWLIB_CXX11
	namespace detail
	{
		template<typename T0, typename... Types>
		struct gcd_return_type
		{
			using return_type = decltype
			(
				  std::declval<T0>()
				% std::declval<
					typename gcd_return_type
					<
						Types...
					>::return_type
				  >()
			);
		};

		template<typename T0, typename T1>
		struct gcd_return_type<T0, T1>
		{
			using return_type = decltype(
				  std::declval<T0>()
				% std::declval<T1>());
		};
	}

	template<typename T0, typename T1>
	inline auto gcd(const T0& a, const T1& b)
		-> typename detail::gcd_return_type<T0, T1>::return_type
	{
		using Type = decltype(a % b);
		return gcd<Type>(
			static_cast<Type>(a), 
			static_cast<Type>(b));
	}

	template<typename T0, typename T1>
	inline auto lcm(const T0& a, const T1& b)
		-> typename detail::gcd_return_type<T0, T1>::return_type
	{
		using Type = decltype(a % b);
		return lcm<Type>(
			static_cast<Type>(a), 
			static_cast<Type>(b));
	}

	/* @brief: 计算 gcd(x0, x1, ..., xn)
	   @complexity: O(nlog m)
	   @param: x0, x1, ..., xn
	   @return: gcd(x0, x1, ..., xn)
	   @return type: decltype(x0 % x1 % ... % xn)
	   @version: since c++11 */
	template<typename T0, typename... Types>
	inline auto gcd(const T0& x0, const Types&... xs)
		-> typename detail::gcd_return_type<T0, Types...>::return_type
	{
		return gcd(x0, gcd(xs...));
	}

	/* @brief: 计算 lcm(x0, x1, ..., xn)
	   @complexity: O(nlog m)
	   @param: x0, x1, ..., xn
	   @return: lcm(x0, x1, ..., xn)
	   @return type: decltype(x0 % x1 % ... % xn)
	   @version: since c++11 */
	template<typename T0, typename... Types>
	inline auto lcm(const T0& x0, const Types&... xs)
		-> typename detail::gcd_return_type<T0, Types...>::return_type
	{
		return lcm(x0, lcm(xs...));
	}
#endif

	/* @brief: 解方程 ax + by = gcd(a, b)
	   @complexity: O(log n)
	   @param: a, b, reference to x, reference to y
	   @return: 返回 gcd(a, b)
	            (x, y) 为 ax + by = gcd(a, b) 的一个根
	   @warning: integer overflow */
	template<typename Type, typename TTemp1, typename TTemp2>
	inline Type exgcd(
		const TTemp1& a,
		const TTemp2& b, 
		Type& x, 
		Type& y)
	{
		if(b == 0)
		{
			x = Type(1);
			y = Type(0);
			return __BWLIB_MOVE(a);
		}

		Type ret = exgcd(b, Type(a % b), x, y);
		Type temp = __BWLIB_MOVE(x);
		x = __BWLIB_MOVE(y);
		y = temp - a / b * x;
		return __BWLIB_MOVE(ret);
	}

#ifdef __BWLIB_WITH_GMPXX
	template<>
	inline mpz_class exgcd(
		const mpz_class& a,
		const mpz_class& b, 
		mpz_class& x, 
		mpz_class& y)
	{
		mpz_class g;
		mpz_gcdext(g.get_mpz_t(), 
			x.get_mpz_t(), y.get_mpz_t(), 
			a.get_mpz_t(), b.get_mpz_t());
		return g;
	}
#endif
}

#endif
